Collective spin systems in dispersive optical cavity QED: Quantum phase transitions 

and entanglement 
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We propose a cavity QED setup which implements a dissipative Lipkin-Meshkov-Ghck model - an 
interacting collective spin system. By varying the external model parameters the system can be made 
to undergo both first-and second-order quantum phase transitions, which are signified by dramatic 
changes in cavity output field properties, such as the probe laser transmission spectrum. The 
steady-state entanglement between pairs of atoms is shown to peak at the critical points and can be 
experimentally determined by suitable measurements on the cavity output field. The entanglement 
dynamics also exhibits pronounced variations in the vicinities of the phase transitions. 

PACS numbers: 42.50.Nn, 42.50.Pq, 03.65. Ud, 73.43.Nq 



I. INTRODUCTION 

The branch of atomic physics associated with ultra- 
cold atoms, ions, and molecules now provides a rich 
and exciting arena for investigations of strongly interact- 
ing, many-body quantum systems. Trapping and cool- 
ing techniques, coherent laser or microwave interactions, 
and applied magnetic fields enable exquisite control of 
both external (motional) and internal (electronic) de- 
grees of freedom of the particles, allowing one to "tai- 
lor" particle-particle interactions and thereby implement 
a broad range of systems that can be described accu- 
rately and transparently by idealized (but nontrivial) 
many-body Hamiltonians. An important example is the 
Hubbard model, realized with ultracold atoms in peri- 
odic optical lattices [Tl 1^, while realizations of other 
novel and significant lattice-spin models have been pro- 
posed, for example, with dipolar molecules in optical lat- 
tices |3] and with chains of trapped atomic ions ^ . The 
common, defining feature of these systems is the possi- 
bility for quantum critical phenomena, i.e., transitions 
between distinct quantum phases, in response to varia- 
tions of an effective field or particle-particle interaction 
strength around some critical value. 

The above-mentioned schemes generally provide many- 
body quantum systems that are subject to short-range 
(e.g., nearest-neighbor) interactions. Another interesting 
and commonly studied class of many-body systems are 
those possessing long-range, or even infinite-range, inter- 
actions, for which theoretical models typically allow exact 
solutions in the thermodynamic limit, or at least enable 
efficient numerical solution for large numbers of particles. 
A standard and classic example is the Lipkin-Meshkov- 
Glick (LMG) model [5], which was originally introduced 
in nuclear physics and is described by a Hamiltonian of 
the form 



i/i 



LMG 



-2hJ,^ 



2A 



(J^+lJy), 



(1) 



where {Jx, Jy, Jz} are collective angular momentum op- 



erators for N spin- 1/2 particles, h and A are parame- 
ters giving the effective magnetic field and spin-spin in- 
teraction strengths, respectively, and 7 G [—1,1] is an 
anisotropy parameter. In this model, each spin interacts 
identically with every other spin and the nature of this 
interaction may be ferromagnetic (A > 0) or antiferro- 
magnetic (A < 0). Significantly, the model exhibits crit- 
ical behavior at zero temperature; in particular, either 
first- or second-order quantum phase transitions may oc- 
cur (depending on the choice of A and 7) as the ratio 
between A and h is varied across a critical value. 

This quantum critical behavior, combined with the rel- 
ative simplicity of the model, has led to renewed theoret- 
ical interest in the LMG model from the point of view of 
studying entanglement properties of many-particle sys- 
tems in relation to quantum phase transitions O |7l [8] . 
Bipartite entanglement measures characterizing entan- 
glement between a pair of spins (e.g., the concurrence) or 
between two blocks of spins (e.g., the entanglement en- 
tropy) are relatively straightforward to compute for the 
LMG model and can display marked critical behavior and 
scaling at quantum critical points [5J ^U\ [TTl [TU [TSJ HH 
[15]. 

Given these interesting and very topical features of the 
LMG model, it follows that the physical realization of 
a system described accurately by such a model would 
provide a valuable test bed for studies of quantum crit- 
ical phenomena and entanglement. However, the ques- 
tion naturally arises as to how realistic such an ideal- 
ized model could be; the assumption of "infinite-range" 
interactions is obviously demanding and implies a very 
specialized system. Hamiltonians of the form ([Tl (with 
7 = 0) have appeared recently in reduced two-mode mod- 
els of atomic Bose-Einstein condensates undergoing tun- 
nelling in double-well potentials or transitions between 
two internal atomic states [TSJ [T7], and in models of a 
few trapped ions interacting with laser fields [THl [T^ . 
but emphasis in these works has been on unitary or adia- 
batic evolution from some initial atomic state to some fi- 
nal, prescribed (entangled) state, while flexibility of these 



systems with respect to parameters of the LMG model 
(i.e., A, N, 7) appears Hmited. 

Another possibihty, furnished by the field of quantum 
optics, and for which long-range atom-atom interactions 
actually occur quite naturally, is cavity quantum elec- 
trodynamics (cavity QED) [20 . Here, one considers en- 
sembles of atoms interacting, through an electronic tran- 
sition, with a common electromagnetic field mode sup- 
ported by an optical resonator. Through this common 
coupling, the field mode can eff'ectively mediate interac- 
tions between atoms located at quite arbitrary and sep- 
arate positions within the mode. So, in particular, the 
concept of an interaction "length" becomes redundant in 
this setting and a collective description of the atoms is 
appropriate. 

In fact, that an ensemble of atoms coupled to a com- 
mon field mode can be viewed as a many-body system 
of interacting spins was highlighted many years ago with 
the prediction of a thermal equilibrium phase transition 
in the celebrated Dicke model of N two-level atoms cou- 
pled to a single quantized field mode [211 l22l [23 l [24 l |25] . 

H ^UJa^a + uJoJz + —^{a^ +a){J+ + J^), (2) 



where a is the annihilation operator for the field mode 
of frequency lo, luq is the atomic transition frequency, 
and A is the atom-field coupling strength (we set h= I). 
In particular, above a certain critical value of the cou- 
pling strength the system enters a so-called "superradi- 
ant" phase [53] • This phase transition persists at zero 
temperature [?7l US], with associated critical behavior 
of both the atom-field and atom-atom quantum entan- 
glement [29l [30] [31]. The critical coupling strength at 
zero temperature is given by Ac — ^^luluq/2, which means 
that A must be comparable to the field and/or atomic 
transition frequencies if the transition regime is to be 
reached. For atomic dipole transitions, this is typically 
not the case and, in fact, if it happened to be so, then 
the model (pi) would be inadequate; in particular, the A^ 
term [omitted from ([2|] of the minimal coupling Hamil- 
tonian should be included and doing so one actually finds 
that no phase transition exists [32J. 

However, a recent proposal for realizing the Dicke 
model quantum phase transition, based on Raman tran- 
sitions between stable atomic ground states in an optical 
cavity QED setting [33], circumvents these issues by (i) 
implementing a system in which the relevant frequency 
and coupling scales are determined by light-induced fre- 
quency shifts and Raman transition rates, and (ii) uti- 
lizing an open-system dynamics (as opposed to a closed, 
Hamiltonian system) with input and output fields (i.e., 
external laser fields and cavity mode losses), thereby re- 
placing a (fragile) thermal equilibrium phase transition 
with a (robust) dynamical, nonequilibrium phase tran- 
sition. Furthermore, as shown in Ref. |33j . the cavity 
output field offers a unique window on the system's be- 
havior and properties, with, for example, fluorescence 
and quadrature-variance measurements providing dra- 



matic signatures of criticality in the system, as well as 
quantitative measures of fluctuations and entanglement. 
These features of the optical cavity QED system, com- 
bined with the observation that, in the dispersive limit 
Lu 3> {wo. A}, the cavity mode may be adiabatically elim- 
inated and the Dicke Hamiltonian reduced to the form 



4A^ 



(3) 



where Jx = |(J+ + J-), motivate us to explore the possi- 
bilities for studying the LMG model in such a setting. In 
particular, by generalizing the conflguration of Ref. [33] 
to two cavity fleld modes and operating in a disper- 
sive regime (amounting to far-off-resonant Raman tran- 
sitions), we find that it is possible to implement atomic 
spin systems that are described by the most general LMG 
model M , and for which the Hamiltonian dynamics may 
still dominate over losses to the output cavity fields, thus 
enabling the clear realization of critical phenomena, in- 
cluding both first- and second-order dynamical quantum 
phase transitions. We find also that the cavity output 
fields can again be used to provide clear and detailed 
probes of properties of the atomic collective-spin system, 
including entanglement, in the critical regime. 

We note that the present work bears some relation to 
studies of optical bist ability and resonance fluorescence 
in cooperative atomic systems, which can also exhibit 
first- and second-order nonequilibrium phase transitions 
(see, for example, [331 [35] [351 [SB [3E] ) • There, however, 
the dynamics explicitly includes (resonant) coherent driv- 
ing of the atomic system by an external laser field (i.e., 
the Hamiltonian describing the system contains a driv- 
ing term linear in J^ or Jy, rather than a direct spin-spin 
interaction term), and relatively little investigation has 
been made of the quantum entanglement associated with 
the critical behavior [39 . 

A more specific outline of our paper is as follows. In 
Sec. [H] we describe the microscopic model of atoms and 
light fields that realizes our effective spin system. In 



Sec. HI we present some background to the LMG col- 
lective spin model and show how to engineer it using the 
general setup presented in Sec. [IT] We conclude Sec. jlllj 
with a brief overview of the methods of analysis to be 
used later in the paper. In Sec. jIVj we describe a more 
specific, potential physical implementation of the system 
we have proposed, based on alkali metal atoms confined 
within a high-finesse ring cavity. In Sec. [V] we focus on 
the 7 = LMG model and focus on the second-order 
transition; we first present a linearized analysis of the 
system in the thermodynamic limit using the Holstein- 
Primakoff representation of spin operators. Using the 
input-output theory of quantum optics we relate the in- 
ternal spin properties to the measurable cavity output 
field and determine the probe transmission spectrum as 
an example. The second part of Sec. [V] is concerned with 
the presence and behavior of atom-atom (or spin-spin) 
entanglement in the system, particularly across the quan- 
tum phase transition (QPT) . We present results for both 



the steady-state entanglement and the entanglement dy- 
namics, using either exact numerical solutions for finite 
system size or analytical solutions in the thermodynamic 
limit. In Sec. |VI| we essentially repeat the analysis of the 
previous section, but focus on a parameter regime where 
a first-order phase transition occurs in the 7 = LMG 
model as the effective magnetic field parameter, h, is var- 
ied. Finally, in Sec. |VII| we conclude and briefly discuss 
possible extensions of the current work. 



II. THEORETICAL MODEL 

We consider a collection of N atoms coupled via elec- 
tric dipole transitions to (at most) four laser fields and 
to a pair of orthogonally polarized optical cavity modes. 
The atomic level and excitation scheme is shown in Fig. IT] 
In particular, the atoms are assumed to possess two sta- 
ble electronic ground states, labelled |0) and |1), at ener- 
gies (S = 1) 0^0 = and wi, respectively. The laser and 
cavity fields combine to drive Raman transitions between 
|0) and |1), via the excited atomic states |r) and \s) (ener- 
gies ojj. and ujs, respectively). Specifically, the laser fields, 
at frequencies ujro, lOsq, tOri, and Wsi, couple to the dipole 
transitions |0) ^^ \r), |0) ^^ |s), |1) ^^ |r), |1) ^^ \s) with 
Rabi frequencies ri^o, f^sOj ^ri, and Q,si, respectively. 
Cavity field a, at frequency w^, couples to the transitions 
|0) <-> |r) and |1) <-^ \s) with coupling strengths g^o and 
gal, respectively, while cavity field 6, at frequency uJi,, 
couples to the transitions |0) <-^ \s) and |1) <-^ \r) with 
coupling strengths gso and g^i, respectively. All of the 
fields will be assumed to be far-off resonance with the 
electric dipole transitions to which they couple, meaning 
that the atomic states \r) and \s) are only virtually ex- 
cited and can be eliminated from the dynamics. Finally, 
at the location of the atoms, the cavity and laser fields 
are taken to be travelling waves copropagating in the x 
direction, with sufficiently broad beam waists so as to 
ensure a homogeneous atom-field coupling. 



A. Adiabatic elimination of atomic excited states 

To facilitate adiabatic elimination of the atomic ex- 
cited states we move to a rotating frame according to 



the unitary transformation U{t) — e 



-iHnt 



with 



Hq = {ujso — uj[)a'a + {ujra — ^'i)b'b 



N 

+ I]Ko|sj)(s,|+w.o|rj)(rj|H-^i|l,)(lj|),(4) 

where uj[ is a frequency close (or possibly equal) to uji. 
Next, as mentioned above, we assume large detunings of 
the light fields from the atomic excited states, i.e., we 
assume that A,. = LUr — cj^o and A^ = lOs — tUso are much 
larger in magnitude than any other rates characterizing 
the system. This allows the atomic excited states to be 




FIG. 1: Atomic level and excitation scheme for the general 
model. 



adiabatically eliminated and also enables us to neglect 
the effects of atomic spontaneous emission. 

Additionally, as depicted in Fig. [l| we assume that 
only four distinct Raman transitions are of significance 
(i.e., resonant or roughly resonant); in particular, in our 
model we retain only those Raman processes that cause 
a change in the electronic state of the atoms (|0) — > |1) 
or |1) — !■ |0)) and also involve transfer of a photon from a 
laser field into a cavity mode or vice versa. All other pos- 
sible Raman processes are assumed to be far-off resonant 
and therefore negligible. Quantitatively, this requires, for 
example, that joj^o ~ i-^al and jw^o ~ (^^^ri +1^1) I are suffi- 
ciently large, with, in particular, |a;,-o — Wa|, \i^rO ~ {^ri + 

Wl)| > \i^a - (t^rl + Wi)|, \uJrO - {uJb + (^l)\- 

Retaining only the four dominant Raman processes 
simplifies the model considerably, and with a choice of 
laser frequencies such that uJsq — u^i = w^o — i-^si — 2w']^ 
we are able to remove all explicit time dependence from 
the Hamiltonian describing our system. Employing the 
collective spin operators. 



1 ^ 



(5) 



j=i 



N N 

J+ = Eii.)(o.i' J--T.\^,)ihl (6) 

and omitting constant energy terms, our effective Hamil- 
tonian for the collective atomic system and cavity modes 



can be written in the form 



H„ 



+ ^{Xaa + Xy) + ^(X,b + xlb^), (7) 



N 



where X, = a^ J+ + Pi J- and the effective parameters 
are given in terms of the microscopic parameters by 



Wo = 



1 /l^.il 
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(8a) 
(8b) 
(8c) 

(8d) 
(8e) 
(8f) 
(8g) 



Note that the (dimensionless) factors {aafi,l3a,b} G 
[—1,1] have been introduced for convenience. Note also 
that for a characteristic level scheme as shown in Fig. [T] 
one might typically expect that gsi = dro and gri = gso, 
so assuming A^ « A^ we would therefore also expect 
that \S-,\ « |<5+,|. 

In summary, the master equation for the reduced den- 
sity operator, pg (i.e., with the atomic excited states elim- 
inated and spontaneous emission neglected), is given by 



Ps 



-■i-[Hg,Pg] + KaD[a]pg + KbD[b]f 



(9) 



where D[A]p = 2ApA^ - A^ Ap - pA^A and Ki is the 
cavity field decay rate. 



B. Adiabatic elimination of the cavity modes 



We now consider the limit ^k| -|- 6f ^ Aa, A(,,ijJo- In 
this limit, the cavity modes are only ever weakly or vir- 
tually excited and may also be adiabatically eliminated 
from the dynamics. Following the standard adiabatic 
elimination procedure |40j . we derive the following mas- 
ter equation for the reduced density operator, p, of the 
collective atomic system alone: 



-^[H, p] + TaD[Xl]p + nD[xl]p, (10) 



with 



H — LOqJz ^-^a^l ^-^bX^ 



N 



^ ^^Y.xl 



where the effective spin-spin interaction strengths and 
collective atomic dissipative rates are (z € {a, 6}) 



A, 

r, 



A?<5, 



AfK, 



(12a) 
(12b) 



Note that both dispersive nonlinear terms [terms propor- 
tional to 5~ and 6^ in Eq. »7l] do not contribute in the 
adiabatic approximation since in this limit we assume a 
vacuum state for both cavity modes. 



C. Cavity output fields and measurement 

Taking a brief step backwards now to the atom-cavity 
Hamiltonian ([7|, and using the input-output theory of 
open quantum optical systems [101 SI], we can derive 
quantum Langevin equations for the cavity mode oper- 
ators; in particular, for the mode b, we have (neglecting 
the term proportional to S^) 



b = — (Kft + i6b)b— iXt 




(13) 



where fein(i) describes the quantum noise input to the 
cavity mode (see Fig. [2| and satisfies the commutation 



relation [bin{t),bl^{t')] = 5{t-t'). Equation (13l illus 



trates the linear relationship between the cavity opera- 
tor and atomic operator Xl . The adiabatic limit of the 
preceding subsection amounts, in the present context, to 
the assumption that X}j{t) varies on a much slower time 
scale than b{t) [and 6i„(i)], so that we can write 



b{t) ~ -i- 



Ab Xl{t) 



Kb + iSb \/N Kb + iSb 



knit). (14) 



The cavity output field is given by 6out(i) — \J2K,b b(t) — 
&in(i), so we in turn obtain a direct relationship between 
the dynamics of the (internal) collective atomic spin and 
the (external) cavity output field. Hence, spin-spin cor- 
relations of the form {XbXb)/N and {XlXb)/N could 
be deduced from correlations of the cavity output field, 
which may be measured, for example, by performing 
broadband homodyne detection on the emitted light |35] . 



N 



III. COLLECTIVE (LMG) SPIN MODELS 

The LMG model, originally introduced in nuclear 
physics to model collective motion in nuclei [5] , describes 
N interacting fermions distributed on two TV-fold degen- 
erate levels (denoted by ±) separated by an energy S. De- 
noting the fermion annihilation operator by c^o-, where 
(11) j G {1, . . . , N} and a € {-I-, — }, the Hamiltonian for this 



system may be written as 



w 
Y 



j,j'-o- 



Et t 



J.-cr^J ,o-- 



(15) 



Introducing the collective spin 



hT>j,a'^c],^Cj,„ and J± = X; 
express the Hamiltonian as 



.t 



J S',±^.J'.=F 



operators, Jz — 
allows us to re- 



H' = 6Jz 



^{Jl + Jl) + ^{J+J- + J-J+).{16) 



This Hamiltonian commutes with J , thus conserving the 
total angular momentum, and with e"''^ , corresponding 
to a parity (spin-flip) symmetry [9] . It is straightforward 
to rewrite this Hamiltonian in terms of J^ and Jy, defined 
via J± = Jx ^ iJy, giving the generalized LMG model. 



H, 



LMG 



= -2hJz 



2A 



{Jl+lJl) 



(17) 



where A = -{V + W)N/2, 7 = (H/ - V)/{V + W) (we 
will only consider 7 g [—1, 1]), and h = —5/2. 

This model is well known for its second-order symme- 
try breaking phase transition in the ferromagnetic regime 
(A > 0) ^3j. For small interaction strength the system is 
in the normal phase, where the ground state is unique and 
polarized in the direction of the magnetic field. As the 
interaction is increased above a critical value, Ac, the sys- 
tem enters the broken phase, where the ground state be- 
comes doubly degenerate and macroscopically displaced 
from its original configuration, thus breaking the parity 
symmetry. For the special case 7=1 the Hamiltonian 
also commutes with J^, thus enabling a direct analytic 
solution. All other cases 7 y^ 1 lie in a separate univer- 
sality class. In the antiferromagnetic regime, A < 0, the 
model exhibits a first-order phase transition as the effec- 
tive magnetic field h crosses /ic = (provided 7 > 0). 

Using the setup described in the previous section we 
can implement the generalized LMG model for any 7 by 
making appropriate choices of aa, f3a,ctb, Pt in the Hamil- 
tonian (111. We now consider three specific cases of gen- 



eral interest. 



A. Conventional 7 = — 1 LMG model 

The 7 = — 1 LMG Hamiltonian may be implemented 
by choosing Ua = Ob — a and Pa — —Pb = P (corre- 
sponding to Xa = aJ+ + /3J_ and Xb = aJ^ — PJ-), 
and setting A^ — —Kb (note that the signs of Aa,b are 
determined by the signs of the detunings 5a,b), so that 



H = -2hJz - 



2A 



{•^x "^y/^ 



(18) 



with h — — wo/2 and A = 2aPKa. This instance of the 
LMG model has been most widely studied for its phase 



transition properties. For the dissipative terms we as- 
sume, for simplicity, that 2Ta = 2Vb = F, so that the full 
master equation reduces to the form 



P = -i[H,p\ 



N 



D[J+]p+^D[J^]p, (19) 



where F+ = Ta^ and F_ = F/?^. The Hamiltonian 
dynamics can be expected to play a dominant role if 

5a,b > Ha,b (which COrreSpOuds to \Ka,b\ > ^a,b)- 



B. Isotropic 7 = 1 LMG model 

The isotropic 7=1 LMG Hamiltonian may be ob- 
tained, for example, by choosing a^ = /?{, = 1 and 
db = Pa = Q{ corresponding to Xa = J+ and Xb = J-), 
and setting A^ = Ab = A, which gives 



H 



-2hJz-'^{j! + j'y), 



(20) 



where h = —ujq/2. The full master equation is 



p = -^[H,p] + ^D[J^]p+^D[J+]p. (21) 



C. Simple 7 = LMG model 

The 7 = LMG Hamiltonian, which will be focus of 
our attention in this paper, may be obtained by choosing 
o-a = Pa = ct (corresponding to Xa = 2aJx), and setting 
5b = Q (so that Ab = 0). This gives 



H 



-2hJz 



2\ 

iV 



r 

•^x^ 



(22) 



where A = 2a^AQ. While the Raman channels involving 
the cavity mode h could be omitted completely, here we 
retain one of them (for reasons to be discussed below), 
with the choice Pb = /3, and Off, = 0, corresponding to 
Xb — PJ-. Hence, the full master equation we consider 
is 



F, 



Tb 



^[H,p] + ^D[2Jx\p+^D[J+] 



N 



N 



(23) 



where the factor /3^ has been absorbed into F^. 

If we now consider the case where |Aa| ^ F^, and 
Fb ^ Fq, then the role played by each cavity mode in 
relation to the atomic system is quite distinct. Specif- 
ically, cavity mode a mediates the collective spin-spin 
interaction required for the Hamiltonian dynamics (with 
coupling strength A^ ~ \^/5a), while cavity mode b ef- 
fectively mediates the collective atomic decay (with rate 
Fb — \1/ Kb). Importantly, we note that Xb = J+ implies 
a quite direct relationship between moments of the cavity 
mode operators b and 5' and moments of the collective 
atomic spin operators J±\ in particular, measurements 
of the output light field from cavity mode h will provide. 



rather directly and transparently, characteristic proper- 
ties of the collective atomic spin. 

In contrast, for the 7 = — 1 model the two cavity modes 
mediate the collective spin-spin interaction on an equal 
footing, i.e., |Aa| = |A(,|, while the operators Xa and 
Xb are linear combinations of J+ and J_, which leads 
to a somewhat less transparent (i.e, arguably less con- 
venient) relationship between correlations of the cavity 
output fields and atomic spin-spin correlations. Partially 
for this reason, we focus in this paper on the 7 = model, 
with a clear distinction between the effective roles of the 
two cavity modes and a potentially better suitability for 
measurements of the collective atomic spin properties. 



D. Methods of analysis 

To analyze the atomic-spin master equations presented 
in the preceding sections, we make use of both numeri- 
cal and analytical techniques. For finite spin j = N/2, 
the master equations can be solved numerically for quite 
large N [43j, owing to the linear scaling of the Hilbert 
space dimension, d, with the system size, i.e., d = iV+ 1. 
In what follows, we will typically present results of nu- 
merical simulations for N < 100. 

For very large system sizes, TV ^ 1, it is possible to 
linearize the quantum fluctuations around the mean spin 
state (i.e., around the "Bloch vector"). First we find this 
mean spin state by calculating the steady-state solutions 
of the semiclassical equations of motion for the compo- 
nents of the Bloch vector. After a suitable rotation (de- 
termined by the mean state) of the spin coordinate sys- 
tem, we use the Holstein-Primakoff (HP) representation 
of angular momentum operators [HJ US] , which enables 
a systematic large- iV expansion of the master equation, 
to which we then apply the limit N ^' 00. While all 
of the results obtained in the linearized regime are exact 
analytical results, in many cases the expressions obtained 
are too lengthy to give any useful information; in these 
cases we simply plot the relevant quantities. 



IV. POTENTIAL EXPERIMENTAL 
IMPLEMENTATION 

For a possible experimental implementation of our 
scheme, we consider, as mentioned previously, an ensem- 
ble of atoms confined inside a high finesse ring cavity 
that supports two travelling-wave modes, a and b. The 
required laser fields, which are assumed to be at frequen- 
cies that are not supported by the resonator, are injected 
through one of the resonator mirrors so as to be coprop- 
agating with the cavity fields through the ensemble. 

If we take ^Li as the atomic species, then the atomic 
level scheme of Fig. [l] can be implemented directly with 
the two ground magnetic substates \F = 1/2, m = ±1/2) 
as |0) and |1), and with a magnetic field applied perpen- 
dicular to the cavity axis to provide a frequency splitting 



2ujb between these two states. The modes a and b would 
be orthogonal, linearly polarized cavity modes, with, in 
particular, mode a polarized along the direction of the 
magnetic field. (Note that if the two modes happen to 
be very different in frequency due, for example, to bire- 
fringence in the cavity mirrors, then the magnetic field 
may not be necessary.) 

Another possibility, illustrated in Fig. |2] might be a 
configuration based on the i^ = 1 ^^ i^' = transition 
of ^^Rb, in which the states |0) and |1) are the ground 
magnetic substates \F = 1,111 = ±1), with frequency 
spHtting 2ujb due to a magnetic field applied along the 
cavity axis. The modes a and b would be orthogonal, 
linearly polarized cavity modes, polarized perpendicular 
to the magnetic field. Note, however, that the modes 
would need to be sufficiently different in frequency (which 
could be imposed, for example, by cavity birefringence) in 
order that the Raman channels involving different modes 
are distinct. 

Alternatively, the modes a and b could be two entirely 
different (linearly polarized) longitudinal modes of the 
resonator, one quasiresonant with the F = 1 ^^ F' — 
transition of the D2 line and the other quasiresonant with 
the F = 1 ^^ F' = 1 transition of the Dl line. 

For specific parameter values, we consider experimen- 
tal systems such as those realized recently in Ref. [46ll47]. 
where cold atoms are held inside a high-finesse opti- 
cal ring cavity. In particular, let us assume a single- 
atom-single-photon dipole coupling strength of 17/(271) ~ 
100 kHz and a cavity field decay rate of Ka/{2Tr) ~ 
25 kHz. For A^ ~ 10^ atoms and a characteristic laser- 
Rabi-frequency-to-detuning ratio of il/A ~ 0.005, we 
have Aa/(27r) ~ 250 kHz {aa = 1)- If we assume a Ra- 
man detuning Sa/{2n) ~ 2.5 MHz > \a/{2TT), Ka/{2Tr), 
we then have, for example, A^ ~ ^a/^a — 2ttx25 kHz and 
Ta ^ Aa{Ka/Sa) ~ 2tt X 0.25 kHz. This illustrates that it 
should be possible to achieve a regime where the (coher- 
ent) Hamiltonian dynamics is dominant over the effective 
dissipation. Note also that, for these parameters, read- 
ily achievable ground state magnetic level shifts {2ujb) 
of tens of MHz would suffice to ensure distinct Raman 
channels. 

The same parameter regime could obviously be cho- 
sen for cavity mode b, but if we consider the 7 = 



model as discussed in Sec. IIIC then we might, for ex- 
ample, assume mode b to be more strongly damped (i.e., 
the two cavity polarizations have different finesses), e.g., 
Kf,/(27r) ~ 250 kHz, and, with smaller Raman transi- 
tion rate Ab/(27r) ~ 25 kHz and detuning Si,/{2tt) ~ 0, 
we would then have Ffc ~ ^l/'^b — 2tt x 2.5 kHz ^ Ta- 
Given these considerations, in the next section, where we 
examine the second-order transition of the 7 = model, 
we will typically employ the set of normalized parame- 
ters {h = 1, Fa = 0.01, Ff, = 0.2}, which give a critical 
coupling strength Ac — 1. 

Finally, we note that the rate for single-atom spon- 
taneous emission (neglected in our model) is estimated 
by T,pfl'^/{4A^) <2ttx 0.04 kHz < Aa,Ta,b, where an 
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(b) 
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as the magnitude of the interaction strength is varied [15] . 
This transition will turn out to be similar to the one re- 
cently studied in the dissipative Dicke model with reso- 
nant atom-cavity interactions (as considered in Ref. [33 J. 
However, it should be noted that in the Dicke model 
the cavity field plays an intrinsic role in the dynamics 
and associated critical behavior, unlike in our present 
model where it has been adiabatically eliminated. Con- 
sequently atom-field entanglement is effectively negligible 
in the present context, while atom-atom entanglement is 
significant and will be the focus of our study. 

In Sec. |V A we consider the spin master equation in 
a linearized regime, appropriate for N ^ 1, and deter- 
mine the transmission spectrum of a weak probe laser. 
Spin-spin entanglement is studied in Sec. |VB| both in 
the thermodynamic limit and for finite N; specifically 
the behavior of the steady-state entanglement, as well 
as entanglement dynamics, is examined in the vicinity of 
the quantum phase transition. 



A. Linearized model 




In this section we study the master equation 



model ( 23 1 in the thermodynamic limit by linearizing 



FIG. 2: (a) Schematic of potential ring cavity system and 
setup for measurement of the output transmission spectrum 
of a weak probe laser field of amplitude £p and frequency i/p. 
(b) Possible atomic level scheme as described in the text. 



atomic exited state linewidth of rsp/(27r) = 6 MHz has 
been assumed. 



the quantum fluctuations around the mean-field state. 
Note that the atom-cavity coupling strengths appearing 
in the effective coupling constants ( [8f| and (8g) scale as 
1/VV, where V is the cavity mode volume. The thermo- 
dynamic limit corresponds to A^ ^ c» and F ^ cxd with 
g == N/V, the atomic density in the cavity, constant. 
Since the thermodynamic limit does not alter the effec- 
tive coupling strengths, which scale as y^, we will hence- 
forth refer to the thermodynamic limit as A^ -^ cx) [28] . 

Firstly, we present the semiclassical analysis which 
determines the mean-field state relevant for A^ 3> 
1. We then expand the angular momentum operators 
around the semiclassical steady state using the Holstein- 
Primakoff representation, thus obtaining a linearized ver- 
sion of the master equation, the eigenvalues of which 
are subsequently analyzed. Finally, we calculate, for the 
linearized model, the transmitted amplitude of a weak 
probe laser through the atom-cavity system as a func- 
tion of the probe frequency, i.e., the probe transmission 
spectrum. This physically measurable quantity probes 
the energy, or eigenvalue, structure of the system and, 
as we will see, provides clear signatures of the dynamical 
quantum phase transition. 



1. Semiclassical equations of motion and steady-state 
solutions 



V. SECOND-ORDER PHASE TRANSITION 

We focus first on the positive field case {h > 0) of 
the 7 = LMG model with ferromagnetic interactions 
(A > 0), for which a second-order phase transition occurs 



The equations of motion for the expectation values 
of the spin components of the Bloch vector, (Jx), (Jy), 
and (Jz), are readily derived from the master equa- 
tion (23 1, but do not form a closed set of equations. 
However, by factorizing all terms (JkJi) — > {Jk){Ji) with 



A), I G {a;, y, z}, which corresponds to neglecting quantum 
fluctuations, we obtain a closed set of equations, which 
we call the semiclassical equations of motion from hereon. 
Introducing the notation X = {Jx)/j, Y = {Jy)/j, 
Z = {Jz)/j, where j = N/2, the semiclassical equations 
of motion are found to be 

X = 2hY-TbZX, (24a) 

Y = -2hX + 2\ZX -TbZY, (24b) 

Z = -2\XY + TbiX^ + Y^), (24c) 

with the constraint X^ + Y'^ + Z"^ = 1 corresponding to 
conservation of angular momentum. 

The steady-state solutions of these equations of motion 
exhibit a bifurcation at a critical coupling strength 



\r = h 



Ah 



(25) 



(note Ac > {h,Tb} for Tb ^ 2h). For A < Ac the stable 
steady-state solutions are 



while for A > Ac they become 

2h 



0, 



^ss — 



-'fss — 



Y — 



A' 



± 

2h 



where 



A = A+.\/A 



A2 - 4/^2 
2AA ' 



^^ 



(26) 

(27a) 

(27b) 
(27c) 

(28) 



The bifurcation at Ac is illustrated in Fig. [3J where, to 
facilitate a comparison between semiclassical and finite- 
N solutions (computed from numerical solution of the 
master equation), we plot the second-order moments 
(J2), {Jy), and {JD (since the finite- iV master equation 
gives {Jx) = {Jy) = for all A). We note that the two ap- 
proaches are already in reasonable agreement for A'^ ~ 50. 



2. Holstein-PrimakojJ representation 

The quantum fluctuations that are neglected in the 
semiclassical analysis can be included in the limit A^ ^ 1 
as a first-order correction. This is achieved by using the 
Holstein-Primakoff (HP) representation of the angular 
momentum operators j44l 145] , which in the present con- 
text takes the form 



N 



Jz 


= y-ctc. 


J+ 


^ ^f^^^ 


J_ 


W^cVl-^ 



(29a) 
(29b) 
(29c) 




FIG. 3: Semiclassical (solid line) and finite-A'^ steady-state 
second-order moments for h = 1, Ta = 0.01, Ft — 0.2, and 
TV = 25 (dotted), 50 (short dashed line), 100 (long dashed 
line). 



where c and c' are bosonic annihilation and creation op- 
erators, respectively, satisfying [c,c'] = 1. In particular, 
if A^ > 1 and (J^) « A^/2, i.e., (c^c) < N/2 (so that the 
Bloch vector points essentially along the z axis), then 
the HP representation of J+ and J_ can be reduced to 
J+ ~ ^/N c and J_ ~ y/N c\ effectively linearizing the 
dynamics. 

In the normal phase (A < Ac), this approach can 
be applied immediately since the steady-state solutions 
ATgg — Yss = 0. However, in the broken phase (A > Ac), 
the steady-state solutions Arss,l^s 7^ 0, i.e., the Bloch 
vector is rotated away from the z axis, and the HP rep- 
resentation is most conveniently applied with respect to 
the new orientation of the Bloch vector. We do this by 
first rewriting the semiclassical steady-state solutions in 
terms of spherical coordinates 9 and (j) as Z^s — cos 9, 
Xss = sin COS 0, and Y^s = sin sin 0, and then apply- 
ing a unitary rotation R = cxp(m • 39) around an axis 
u = (— sin0, cosi/), 0), so that the transformed operators 
J'l = R} JiR describe quantum fluctuations around the 
semiclassical steady state. The HP representation (29a)- 



( 29c I and subsequent large- A'^ expansion is then applied 



to the operators {J/}. 

The master equation obtained in this way may be writ- 
ten, for both phases, in the general form (omitting con- 
stant energy terms in the Hamiltonian) 



-Km 
-iV 






-i[Hu,„ p] + T+,kD[cl]p + T^MD[ck]p 



-,fe 



2ckpck + 2clpcl~{cl + icl)^p} 
-2ckpck + 2clpcl - {-4 + (4)2, p} 



(30) 



with 



Hlin = Ai l.ctcfc + A2,k 



J^2 



.t^2 



+ *^3,fc (4) 



4 + icl) 



Cfe 



(31) 



where k G {<,>} and c< (c>) denotes the bosonic op- 
erator for the normal (broken) phase. The coefficients in 
the normal phase are given by 



(32a) 
(32b) 
(32c) 
(32d) 
(32e) 
(32f) 
(32g) 



Ai,< 


= 2/i-A, 


A2,< 


= -A/2, 


^3,< 


= 0, 


r+,< 


= r„. 


r-,< 


= Ta+Tb 


n,< 


- r,. 


ri,< 


- 0, 



while in the broken phase they are given by 
1 



A2.> 
^3,> 

r±.> 



'4h^ - 3Tl + 4AA) 



2A ' 

1 
4AA 

Tb 
4AA 



{Tl-Ah')Jx^-Tl-ihTl]i33h) 



Ah^ + Tl+ihJX^-Tl] , (33c) 



4A2 



r„ 



{Ah 



2 rg; 



A2 - rg + 4hTl 



2A2A 

r 

4TA2 



+ 7^^A^^(4/.A,-2AA), 



p2 

^ b /a2 



4AA2 



(A"-4/i"). 



(33a) 



(t2/i + A)2, (33d) 



(33e) 



(33f) 



u = Mu, where u = ((cfc), (c|,))"^ and M is a 2 x 2 ma- 
trix. The real and imaginary parts of these eigenvalues 
are plotted in Fig. |4]for our characteristic set of numerical 
parameters. We note that except for the region near zero 
coupling strength the eigenvalues exhibit very similar be- 
havior to that found in the dissipative Dicke model [33] . 
In the normal phase (A < Ac) the eigenvalues of M are 
given by 



/i± = — Fb ± 2i^Jh{h — A) , 



(34) 



the imaginary parts of which go to zero at the point A' = 
/i < Ac, with a characteristic scaling of ^/X' — A. For 
A' < A < Ac the eigenvalues are real and distinct, with 
one going to zero at Ac (i.e., critical slowing down) and 
the other to — 2F;,. 

In the broken phase (A > Ac) the eigenvalues of M are 
given by 



M± 



2Tbh 

'~A 



± J2{2h^ 



n 



AA) 



(35) 



In the region A > A", where A" = (Fg + 2h'^)/y/4J}X, 
the eigenvalues are complex conjugate pairs with a real 
part that diminishes for A ^ Ac like —Tbh/X. Provided 
then A" 



Tb < \/2h\/l + \/5, then A" > Ac and the imaginary 
parts vanish as A approaches A" from above with the 
scaling \/X — A" (which can be shown using a first-order 
Taylor series expansion about A = A"). The imaginary 
parts are zero in the interval Ac < A < A", while the 
real parts again approach and — 2Ffc, respectively, as 
A ^ Ac. 

If Tb > \/2h\J\ + \/5 then A" < Ac, and the eigen- 
values are complex conjugate pairs immediately above 
the critical point. In this situation, the dissipation is 
stronger than the Hamiltonian dynamics; this is also an 
interesting regime, but not one that we will consider in 
the present paper. 



Note that the Hamiltonian (31 1 does not contain any 



terms linear in c^ and c\. , which is a consequence of the 

applied rotation, and also means that {cf;)ss — (c|,)ss = 0. 
While the coefficients for the broken phase are rather 
complicated, they do simplify considerably in the limit 
of very large A; in particular, for A 3> h^Ta^b one finds 
Ai,> ~ 2A, A2,> ~ 0, and A3,> ~ 0, while F±,> ~ Fb/4, 
F^ y ~ — Fb/4, and Fl > ~ 0. The master equation then 
corresponds to that of a simple quantized harmonic oscil- 
lator coupled to a somewhat unconventional (squeezed- 
type) reservoir [15] . 



3. Eigenvalue analysis 

It is interesting to examine the eigenvalues associated 
with the linear set of equations of motion for the first- 
order moments (cfe),(cl), which may be expressed as 



4-. Probe transmission spectrum 

A standard way to examine the structure and dynamics 
of an atomic system is to measure the transmission of a 
(weak) probe laser field through the medium as a function 
of the probe frequency. This amounts simply to detecting 
the frequency response of the system to an applied field 
or "force". A schematic diagram illustrating the setup 
for such a measurement in the present context is shown 
in Fig. I (a). 

For our theoretical investigation of the transmission 
spectrum we retain the two cavity modes in our model 
and make use of the input-output theory of open quan- 
tum systems [40l |4T]. In particular, our starting point 
is the atom-cavity Hamiltonian (pi and we again con- 
sider the limit N :$> I, so that we can perform a lin- 
earization. To do this, we follow our previous working 
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FIG. 4: Eigenvalues of the lin ear ized equations of motion, 
fi±, as given by Eqs. (34 1 and (351, for h = 1 and Ft, — 0.2. 



The right-hand column gives a magnified view of the region 
around Ac — 1.01. 



and determine the stable semiclassical steady-state am- 
plitudes of the atom-cavity system from the semiclas- 
sical (i.e., factorized) equations of motion for the mo- 
ments {{a),{b),{Jx),{Jy),{Jz)}- Note that the steady- 
state cavity mode amplitudes in this approach can be 
expressed in terms of the atomic amplitudes as (neglect- 
ing terms proportional to S^^^ and setting Sb = 0) 



N 



-2iXa ^ {b)s: 



A. 



Kb 



(Fs, - iX,., 



(36) 



Using the HP representation of the atomic spin operators 
and linearizing about the semiclassical steady states as 
before leads to the following Hamiltonian for the normal 
and broken phases. 



H, 



g>lin 



Scclck + Sattlttk + 6bb\hk 

+{BiCk 



B2cl)bk + {Blcl 



B*2Ck)bl (37) 



where k e {<,>}, ak and bk denote the annihilation op- 
erators for the intracavity modes in the normal and bro- 
ken phases, and the coefficients {Sc,A,Bi,B2} are given 
in Appendix \K\ 

Employing the quantum Langevin equations of the 
input-output theory of open quantum systems we can 
analytically solve for any cavity output correlations and 
spectra of interest [3S]. Here, however, we focus sim- 
ply on the amplitude of a probe laser field transmitted 
through the system and into the output field, as depicted 
in Fig. [2] (a) . We consider only the case in which a probe 
laser of frequency Vp (in the rotating frame) and ampli- 
tude £p drives cavity field mode b. 

The analytical expression for the amplitude of the 
transmitted probe, Ap(i>), is rather complicated, but if 



we restrict ourselves to a frequency range where |z^| ^ 
5a, Kb (also with Ka <C Sa), then for A < Ac the transmit- 
ted probe intensity is well approximated by 



T,(:.) = |A»f 



iVh 



A^h{h^\) 



u-2y/h{h- \) 



iTh 



^-y/h^^ 



v + 2^Jh{h~ A) +iTb 



(38) 



where we have normalized the intensity such that it takes 
a maximum value of unity for an empty cavity. This form 
for Tp{v) highlights the presence of atomic resonances at 
the frequencies v = 1tii{^,±), superimposed on a broad 
background corresponding to the bare cavity mode reso- 
nance. Note that this is in contrast to the findings in the 
dissipative Dicke model [331 where the probe laser trans- 
mission spectrum exhibits strongly coupled atom-cavity 
resonances. 

In Fig. ^ we plot the transmission spectrum [computed 
from the full theory - note that the approximate expres- 
sion ( |38| is in good agreement for the parameters chosen] 
for a series of values of A on either side or the transition. 
Note that we have chosen Vb — 0.05 here in order to 
highlight the main features of the spectrum more clearly. 
For A <C Ac, we observe, at i^ ~ 2/i, a single sharp dip of 
width 2Tb in the envelope of cavity mode resonance, cor- 
responding to a cavity-mediated collective atomic emis- 
sion resonance {5^ — 2/i) ; at this A, spin-spin interactions 
mediated by cavity mode a are small and have little effect 
on the spectrum. 

As the interaction strength A is increased, this dip 
moves to smaller frequencies and reduces in depth (even- 
tually inverting) , while a peak emerges at the correspond- 
ing negative frequency. The positions and widths of these 
features reflect the real and imaginary parts, respectively, 
of the eigenvalue structure of the system, while their "in- 
tensities" also relate to the populations of the energy lev- 
els. At \ — h the two peaks merge into a single peak cen- 
tered at z^ = 0, with a height Tp{Q) ~ [h/Tbf. Then, as 
A ^ Ac, this peak diverges (corresponding to eigenvalue 
/i_ ^ 0) in a pronounced signature of the phase transi- 
tion. A similar divergence in the probe laser transmission 
spectrum is found in the dissipative Dicke model |33j . 

Just above the critical point, two peaks reappear in the 
spectrum and move apart with increasing A, as shown 
in Fig. Is] The negative frequency peak diminishes in 
strength, while the peak at positive frequency inverts to 
a dip, which narrows and moves to increasingly larger fre- 
quencies. In fact, for A 3> 1, its position is approximated 
by 2A and its width by 2Tbh/X. 
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FIG. 5: Transmission spectrum in the linearized regime, for 
A = 0.3,0.93,0.992, 1.000625(= Ac), 1.005, 1.05, 1.5, with mi- 
croscopic parameters Ka ~ 0.3, Sa ~ 15, and Xt — 0.87, 
Kb — 15, giving Ft = 0.05. We set h — 1 as usual. Note 
that Xa is chosen to give the indicated A for the given choice 
of Ka and da , viz. Eq. ( |12a| and recalling that A = 2Aa , while 
Fa varies according to Eq. (|12b|), with Fa — 0.01 when A = Ac. 



B. Entanglement 

1. Entanglement criteria 

Recently a criterion for bipartite entanglement in col- 
lective spin systems was derived ^50j , and the connection 
to spin squeezing established. For the case of symmetric 
states the criterion is both necessary and sufficient, and 
reads 



a 



N 



(AJ; 



N^ 



W > 0, 



(39) 



where J^p — sm(ip)Jx + cos{ip)Jy. In this work, we will 
use the magnitude of C^ as a quantitative measure of the 
entanglement in the system. Note that for finite N, and 
also in the linearized analysis, we have (J^p) = [since 
there are no linear driving terms in the effective Hamilto- 
_ ^ l-(4/7V)(j2). Note 

also that'Ci^^o = Cy, which was shown to be equivalent 



nians (22 1 and (31 1], and thus C^ 

3h wf 
C [51], in nondissipative LMG mod 



to the concurrence, 
els 13 . 

We also compute the rescaled concurrence, Cr — 
{N — 1)C, which is the relevant (nonvanishing) quantity 
to study for infinitely coordinated collective spin systems 
in the thermodynamic limit [121112] ■ It is possible to show 
that for the system considered here, Cr may be written 
as [52l 



Cr = 



2max{0,Ci} if £■ < F 
2max{0,C2} if £■ > F 



where 
C2 - 



(Jl 



{Jl 



(J. 



N N 2 ' 

4 N 
V[(7V(iV - 2) + 4( J|)]2 - [4(JV -JKhW 
4iV 



(40) 



(41) 



and 



E = 



F 



N 
~2 



N 



(42) 



V[(iV(^-2)+4(J|)]2-[4(7V-l)(J,)]^ 



4:N 



iJi 



N 



(43) 



As pointed out in Sec. |IIC| the spin variances required 
to compute the entanglement measures described above 
can in principle be determined from appropriate mea- 
surements performed on the cavity output field. 



2. Steady-state entanglement 

For finite N we numerically solve the master equation 
for the steady-state density matrix and then compute the 
operator averages required to determine C,^ and Cr. In 
Fig. |6] we plot max{0, C^} as a function of A and <y9 for 
N — 100. We see that below the critical point, A < Ac, 
entanglement is present for a broad range of angles (p. 
However, as the critical point is approached the range 
of angles ip which gives nonzero entanglement, dp > 0, 
becomes increasingly narrow. Once above the transition, 
A > Ac, the region of finite C^ continues to narrow until 
it eventually disappears altogether. 

To help interpret the behavior of C^ , we make use of a 
phase space representation of the atomic state that em- 
ploys the spin coherent states, which are defined by [52] 



\v)^a+\vn 




J+ni 



\j,m)j, (44) 



e'"^ tan ^ , with 9 and corresponding to spher- 
and 



where rj 

ical coordinates, and \j, m) are the Dicke states with 
m e [-j, -j + l,...,j - 1, j] (for our system, j = A^/2). 
Using these states we can define the spin Q-function, 



3(^) = {'n\p\-n)- 



(45) 



Fig. It] displays Qs{v) on the surface of the Bloch sphere 
for A^ = 50 and for a series of interaction strengths A. 
Below the critical point, Qs{r]) is single-peaked and cen- 
tered around the top of the Bloch sphere (9 = 0), with 
little obvious angular dependence. Correspondingly, the 
entanglement measure C^p is finite over a rather broad 
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FIG. 6: Entanglement measure max{0,C^} for A^ — 100, 
/i = 1, Ta = 0.01, and Ft, = 0.2. 



comes increasingly elongated along a direction close to 
the X axis, until, at the transition, it splits into two peaks 
located approximately at the two semiclassical steady- 
state amplitudes (27b I and (27c I. These peaks con- 



tinue to move apart in phase space as the interaction 
strength is increased further; eventually both peaks will 
lie in the equatorial plane corresponding to 6 — 7r/2 and 
4> — 0, TT. Correspondingly, the range of ip over which 
C^ remains finite becomes increasingly narrow and is fo- 
cussed around an axis perpendicular to that along which 
the two peaks lie. This narrowing of the "width" of C^p 
can be explained by noting that, since ( J^) = 0, we have 
C^ — I — (4/iV)([sin(iy9) Jj; -I- cos{ip) Jy]'^) . For increas- 
ing interaction strength A > Ac, (J^) becomes of order 
j^ = N^/A (see Fig. ^, and so the optimal choice of ip 
becomes more critical. In fact, one can show for A > Ac 
that the range of (p over which max{0, C^} > scales as 

1/Vn. 

Next, we consider the rescaled concurrence, Cr, as a 
function of the interaction strength A. In fact, one finds 
that 



range of (p, with a maximum close to 93 = (i.e., near 
Cy). Note that this slight shift of the optimum away 
from (/? = is a consequence of the dissipation (F;,) in 
the system. 



Cr = max Cif 



(46) 




FIG. 7; (Color online) Steady-state spin Q-function, Qsijf), 
on the Bloch sphere for (a) A = 0.5, (b) A = 1.01, (c) A = 1.1, 
and (d) A = 2, with A^ = 50, ft = 1, Fa = 0.01, and Vb = 0.2. 
Note that dark blue corresponds to the minimum value of 
zero of Qs (rf) while dark red indicates the maximum value of 
Qs{ri). 

As A increases towards the critical point, Qs{'n) be- 



i.e., Cr is simply the optimal value of the quantity dp 
just considered. In Fig.[8]we plot Cr versus A and observe 
that the entanglement reaches a maximum for A close to 
Ac [at finite N the critical point is slightly shifted from Ac 
as given in Eq. (25 1]. This peaking of the entanglement 
at the quantum phase transition has been conjectured 
and demonstrated theoretically for the equivalent closed 
(nondissipative) spin models [HI [HJ [T5]. Our results 
confirm that this behavior can persist in steady state in 
our nonequilibrium, open-system version of these models, 
and can in principle be measured within our proposed 
setup. 




FIG. 8: Rescaled concurrence Cr for A'' = 100 (dashed line) 
and in the thermodynamic limit (solid line) with h = l,Va = 
0.01, and F;, = 0.2. 
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In the linearized treatment {N ^ 1) of the HP rep- 
resentation, we can write J^p « {\'N/2)Xip, where Xv — 
i{-Cke''^ + cle-"^). Noting that (cfc) = 0, the entangle- 
ment measure C<^ can be expressed as 



C, 



HP 



1 - ixl) 



[e'^'^{cl) + e-'^^{{cir)-2{clc,)), (47) 



while the rescaled concurrence is Cr is given by 



Q 



HP f 2max{0,Ci"P} if E^^ < F 
^ ~ { 2max{0,C2"P} if E^^ > E 



HP 
HP 



where 



cr - \{ci)\-{cic,), 



Cr = {cW) ^j{{clc,Y) (4.c,> 



and 



i?"P = 2(4c,), 



(48) 
(49) 



(50) 



;.HP 



\J{{c{ck?)-{cick) + \{cl)\. (51) 



Using the linearized master equation ( 30 1 , we can derive 



a closed set of equations for the second-order moments 
(CfeCfc), (c^.), and ((cj,)^), from which we may determine 
the steady-state solutions analytically. Note that the 
fourth-order moment appearing in F^^ can be expressed 
in terms of second-order moments, since the states we are 
dealing with in this linearized approximation are neces- 
sarily Gaussian. 



■< 1 




FIG. 9: Entanglement measure max{0,C^} in the thermo- 
dynamic limit for h = 1, Ta = 0.01, and F;, = 0.2. Above 
the critical point, the system is linearized about only one 
of the two semiclassical steady-state amplitudes; hence the 
less-sensitive dependence of max{0,C^} on ip for A > Ac as 
compared with the finite- A'' results. 



In Fig. l9] we plot C^ as a function of ip and A, as 
determinea from the linearized HP representation. The 
behavior below the critical point (A < Ac) is very sim- 
ilar to the finite-A'' case. However, the behavior above 
the critical point (A > Ac) is very different. Here, the 
sensitivity of C^p to Lp, for A > Ac, is much less critical 
because the linearized model describes only the fluctu- 
ations around one of the two semiclassical steady-state 
amplitudes (i.e., around one of the two lobes appearing 
in the spin Q-function for A > Ac). 

Note that we can obtain plots of max{0, C^} similar 
to Fig. [6] for the region A > Ac, but determined from 
the linearized HP model (with a finite value of A^), by 
making a rotation back to the original coordinate system 
and then setting, by hand, (xv?) = 0; to mimic an equal, 
incoherent mixture of the states associated with the two 
semiclassical amplitudes. 

Finally, returning to Fig. [8] we have plotted C-r as a 
function of A, computed from the HP model in the ther- 
modynamic limit. Again, Cr corresponds to the value 
of Ctp optimized over if, and, since the optimal (f corre- 
sponds to an axis perpendicular to the (above-transition) 
splitting of the semiclassical amplitudes, we expect, and 
indeed find, good agreement with the finite-A^ results 
over the full range of A. 

If we make the simplifying assumption that Fq ~ 0, 
then, for A < Ac one can show that E^^ < E^^ and 



a 



HP 



A(V4/i(Ac-A) + A2- A) 

4/i(Ac - A) 
1 1 /i(Ac - A) 



2 2 A? 



(52) 
for Ac-A<A. (53) 



This shows reasonable agreement with the plot, but 
reaches a maximum value of 0.5 at the critical point. 



3. Entanglement dynamics 

We now consider the dynamics of the entanglement; 
starting from an initially unentangled state, we exam- 
ine the time evolution of the state and its entanglement 
as quantified by CR(t). The initial state is taken as the 
A = ground state, i.e., the state with all atomic spins 
pointing up (which is a convenient state to prepare in an 
experiment). As in the previous section, we compute the 
entanglement both numerically for finite A^ and analyti- 
cally for A^ 3> 1 in the linearized approximation, 
we plot CR(i) 



In 



A^ 



Fig. 
100. 
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versus A and time t for 
At long times we recover the results of the 
previous section, but at short times the behavior as a 
function of A is quite different; the entanglement rises 
to a high value and remains at that value for increasing 
interaction strength A. This behavior can be attributed 
to the Hamiltonian dynamics, which dominate the dis- 
sipation at short times and can create highly entangled 
states. The potential of such Hamiltonian dynamics for 
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generating such highly entangled states has been pro- 
posed previously, for example, in Refs. ^7\ [18]. Note, 
however, that the presence of the term —2hJz in our 
system Hamiltonian tends to make the generated states 
more complicated and less straightforward to interpret. 
Although the focus of this paper is on the quantum phase 
transition, it is clear that with a slight modification the 
scheme also has interesting potential for the controlled 
generation of specific, highly entangled multiatom states 
(e.g., Greenberger-Horne-Zeilinger states) . In connection 
with this, an important aspect of our implementation 
should be highlighted here: because both the effective 
interaction and dissipation of the spins is controlled by 
the optical laser fields, we can in principle "freeze" the 
state of the atomic system at any time by simply turning 
these fields off. 



-•^ 10 





FIG. 10: Rescaled concurrence Cuit) as a function of A and 
t for iV = 100, /i = 1, Fa = 0.01, and Ft = 0.2. 



FIG. 11: Rescaled concurrence Cnit) as a function of A and t 
in the thermodynamic limit, with /i = 1, Fa = 0.01, and F;, = 
0.2. Above the critical point the dynamics is linearized around 
one of the two possible semiclassical steady-state amplitudes. 



icd, one occurring at positive h (the transition discussed 
in the previous section) and the other, equivalent tran- 
sition occurring at negative h. However, in this section 
we show that with the addition of dissipation this model 
actually exhibits a first-order phase transition near /i ~ 
(note that in the absence of dissipation no such transi- 
tion exists). As in the previous section, we begin with 
a study of the linearized spin master equation, includ- 
ing an eigenvalue analysis and calculation of the probe 
transmission spectrum, after which we focus again on 
the entanglement properties of the system. For numeri- 
cal calculations we will typically employ the set of nor- 
malized parameters {A = 1, Fa = 0.01, F(, = 0.2}, which 
correspond to a critical effective field strength /ic — 0. 



In the linearized regime, TV ^ 1, we solve the equations 
of motion for the second-order moments, (c|.Cfc), (c^), 
and ((cj.)^), with the initial conditions (c[.(0)cfe(0)) — 0, 

(4(0)^) = (cfe(O)^) = 0. The results for CR(i) are shown 
in Fig. |ll[ Below the critical point the behavior is similar 
to that observed for finite N. However, above the critical 
point, where the dynamics is linearized about only one 
of the two allowed semiclassical steady-state amplitudes, 
the rescaled concurrence is, as expected, quite different, 
owing to the more limited range of entangled states that 
the linearized (Gaussian) theory can accommodate. 



A. Linearized model 

As before, we consider the thermodynamic limit and 
linearize the master equation ( 23 1 about the mean-field 



state. To do so, we first find the semiclassical steady- 
state solutions and then expand the angular momentum 
operators around these mean-field solutions using the 
Holstein-Primakoff representation. 



Semiclassical steady-state solutions 



VI. FIRST-ORDER PHASE TRANSITION 

We now turn to the case of a fixed, positive interaction 
strength ( A > 0) of the 7 = model (Sec. Ill C I with vari- 
able effective field h. In the absence of dissipation this 
model exhibits two second-order transitions as h is var- 



From the (factorized) semiclassical equations of motion 
for X, Y, and Z, Eqs. ( |24a l-(24cl we again obtain the 
stable steady-state solutions. These exhibit discontinu- 
ities at the critical field strengths 



he. = 



X-JX^-Tt 



(54) 



and h = 0. For h < the stable steady-state so- 
lutions are given by Eq. ( 26 1 , while for he < h < 



(A + ^A ^ — r^)/2 [5i] the stable steady states are given 
by Eqs. (|27a|)-p7cl). While outside the region < h < he 



the stable steady states are unique, inside the region 
< h < he both steady-state solutions (pGJ) and (27a)- 
(27c) are in fact stable. However, for the characteris- 



tic parameters we consider here this region is very small 
{he ~ 0.01). Moreover, we have verified (from a lin- 
earized analysis) that the steady-state solution (26) is 
more stable in the region < h < he and thus we will 
only consider this solution in that region. Note that for 
larger values of the dissipation, Ff,, this region becomes 
more pronounced (in this case all stable steady states 
should be considered [53; ) , but this is beyond the regime 
we wish to consider here. 

The relevant stable steady-state solutions Zgs , -Ygs and 
Yss are plotted in Fig. [12] together with results from nu- 
merical solutions of the master equation for a range of 
values of N up to 100 (at which agreement between the 
two approaches is already quite good). The discontinu- 
ous jump of Zss at /ic — signifies the first-order phase 
transition. Note that for the case A < the same first 
order transition occurs except that it is shifted to —he 
(i.e., in Fig. 12 all curves are flipped about h = 0). 
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FIG. 12: Semiclassical (solid line) and finite-A^ steady-state 
solutions for A = 1, r„ = 0.01, Fb = 0.2, and iV = 25 (dotted), 
50 (short dashed line), 100 (long dashed line). Note the inset 
in the bottom panel is a magnified plot of the semiclassical 
solution of {J'y) I f" ■ 
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fc =<, while for h > he (broken phase) the linearized 
master equation is also identical to Eq. ( 30 1 but with 



3. Eigenvalue analysis 

The eigenvalues of the linearized system, i.e., of the 
matrix M, where u = Mu and u= ((c), (c^))'^ , for h < he 



Eq. (35 



are given by Eq. (34) while for h > he they are given by 
In Fig. |13| the real and imaginary parts of the 
lues are plotted for our characteristic set of pa- 
in the normal phase {h < he) we see that the 

In 



eigenvar 
rameters. 

imaginary parts go to zero at the point h' = < he 
the region h' < h < he both eigenvalues are real and dis- 
tinct, with one going to zero at he and the other going 
to — 2Ff,. This behavior is the same as that found for the 
second-order transition of the earlier section. However, 
immediately above the transition, h > he, the eigenvalues 
become complex conjugate pairs with a nonzero real part 
that diminishes for h^ he like —T^h/X. This discontin- 
uous jump of the eigenvalues is an additional signature 
of the first-order phase transition. 



c 




0.02 



0.02 



FIG. 13: Eigenvalues of the linearized equations of motion, 
fji,±, as given by Eqs. (341 and (351, for A = 1, and Fj, = 0.2. 
The right-hand column gives a magnified view of the region 
around h^ — 0.0101. 



4- Transmission spectrum 



2. Holstein-PrimakojJ representation 

Here, we again include the quantum fluctuations for 
A'' ^ 1 as a flrst-order correction by linearizing the spin 
operators around the semiclassical steady state via the 
HP representation. For h < he (normal phase) the lin- 
earized master equation is identical to Eq. (|30| with 



We determine the probe transmission spectrum in the 
linearized regime following exactly the same calculations 
as outlined in Sec. IV A 41 The linearized Hamiltonian 
describing the full atom-cavity system is easily obtained; 
for h < he it is given by Eq. (37) with 6 = and (f> = 0, 
while for h > he it is also given by Eq. (37), but with 
6 and (t> given according to the semiclassical solutions 



Eqs. (27a)-(27cl as explained in Sec. VA2 
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FIG. 14: Transmission spectra in the linearized regime, for 
h = -0.6,-0.1,-0.01,6.25 x 10"''(= /ic), 0.05, 0.3, with mi- 
croscopic parameters Ka = 0.3, Sa = 15, \a = 2.7, Aj, — 0.87, 
and Hit = 15, giving A = 1, Fa = 0.01 and Fj, = 0.05. 



critical point is approached from below, significant entan- 
glement persists, but for a somewhat narrower range of 
angles (p. However, immediately above the critical point 
the entanglement drops suddenly to zero for all values of 




Restricting ourselves again to a frequency range where 
\v\ <ti \Sa\,Ki, (also with Ka ^ \Sa\), then as previously an 
approximate expression for the transmitted probe inten- 
sity can be derived in the normal phase for h < h^ and 



takes exactly the same form as Eq. ( 38 ) 



In Fig. UM we plot the transmission spectrum for a se- 
ries of values of h across the critical point he. In the 
normal phase {h < he), we observe the same behavior as 
in the normal phase of the system in the previous sec- 
tion (A < Ac), except that the orientations of the peaks 
and dips have inverted in accordance with the change of 
sign of the field {h < 0). The central peak diverges as 
the critical point is approached from below in the normal 
phase, again signifying the phenomenon of critical slow- 
ing down in the vicinity of the phase transition. However, 
immediately above the critical point, the spectrum splits 
discontinuously into two sharp peaks of width ^ Thh/X, 
located at frequencies v ~ ±2A. This jump from a single 
divergent peak at i^ = to a two-peaked spectrum of- 
fers a pronounced, observable signature of the first-order 
transition. 



B. Entanglement 

1. Steady-state entanglement 

We compute, as before, the entanglement measures C^ 
and Cr, both numerically for finite N and analytically 



for A^ ^ 1 in the linearized regime. Fig. 15 shows a plot 
of C^ as a function of h and (p ioi N — 100. We see that, 
well away from the critical point, substantial entangle- 
ment is present over a broad range of angles if. As the 



FIG. 15: Entanglement measure max{0,C^} as a function of 
h and ^p ioT N = 100, A = 1, F„ = 0.01, and Ft = 0.2. 

To help understand these results we again utilize the 
atomic coherent state representation and study the spin 
Q- function Qsiv)- In Fig- 16 we plot Qsii]) on the 
Bloch sphere for a series of values of h in the vicinity of 
the first-order transition. Well below the critical point, 
in the normal phase, Qs{il) is a single peaked function 
with little angular dependence. Correspondingly, C^ is 
nonzero over a broad range of (p, with a maximum close 
to (fi = 7t/2 (i.e., near Cx)- Again, note that this slight 
shift of the optimum away from ip = tt/2 is a consequence 
of the dissipation (Ft,) in the system. 

As h increases towards the critical point, Qs{ti) be- 
comes increasingly stretched along the y axis. As the 
critical point is traversed Qsiv) rapidly rotates around 
from the y axis towards the x axis, and splits into the fa- 
miliar two-lobed structure associated with the two semi- 
classical steady-state amplitudes of the broken phase. At 
the same time as the critical point is approached, the 
range of p over which C^ remains finite narrows and im- 
mediately above the critical point it drops abruptly to 
zero for all choices of p. This behavior is akin to the be- 
havior we observed for large interaction strength in the 
regime of the previous section, where (Jj) becomes of or- 
der j^ = N'^/4 (see Fig. 12 1, which severely restricts the 



range of p for which Cip > 0. Note that at larger values of 
h than displayed in Fig. 1151 the entanglement, C^, once 
again becomes nonzero (centered around (/9 « 0) coincid- 
ing with the broken phase behavior of the second-order 
transition discussed in the previous section. 

In Fig. [17] we plot the rescaled concurrence Cr as a 
function of the effective field strength h and again find 
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FIG. 16: (Color online) Steady-state spin Q-function, (5s('7), on the Bloch sphere for (a) h = —0.5, (b) h — —0.01, (c) 
ft = 2.5 X 10"^ (d) /i = 5 X 10~^ (e) h = 0.015, and (f) h = 0.15, with N = 50, X ^ 1, Fa ^ 0.01, and Ft, = 0.2. Note that 
dark blue corresponds to the minimum value of zero of Qsi'n) while dark red indicates the maximum value of Qsiff)- 



that close to the critical point, h^, the entanglement 
reaches its peak value. Although the equivalent closed 
system would not feature a maximum in the entangle- 
ment near h^ (due to the complete absence of a phase 
transition), this result is in agreement with a conjecture 
concerning entanglement in open systems at quantum 
critical points [55] . 
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FIG. 17: Rescaled concurrence Cr versus h for A'^ = 100 
(dashed line) and in the thermodynamic limit (solid line) with 
A = 1, Fa = 0.01, and Ft = 0.2. 

In the linearized treatment {N ^ 1) we obtain very 
similar plots of C^ to those of finite N (Fig. 15 1 and for 



Cr the result is shown in Fig. |17[ In the limit where we 
consider Fq ~ we can again obtain an approximate ex- 
pression for the rescaled concurrence (for h < he) given, 
in this instance, by 



a 



HP A(V(fe-A/2)(fc- fec) + A2-A) 

- h,) 

he- h't: X. 



A{h-A/2){h 

1 Ihe-h 

2 " 2 A 



for 



(55) 



This again has a maximum value of 0.5 at the critical 
point, and, for large \h\, drops off like l/|/i|, in reasonable 
agreement with the plots. 



2. Entanglement dynamics 

Finally, in Fig. [18] for N = 100 we illustrate the time- 
dependent behavior of the rescaled concurrence, CR,(i), 
for varying h, given an initial (unentangled) state with 
all spins up. Once again, we observe an interesting os- 
cillatory behavior of C^(t), with, in particular, highly 
entangled states generated by the Hamiltonian dynamics 
at short times (for almost all values of h), before dissi- 
pation has had time to play a significant role. For the 
linearized regime {N ^ 1) a similar plot of CR(t) can be 
obtained which agrees well with the finite N result for 
h < he but shows zero entanglement for almost all values 
of h > he because of the restricted linearization around 
only one of the two permitted semiclassical steady-state 
amplitudes. 
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amining a system of multiple (separately addressable) 
atomic pseudospins all coupled to the same quantized 
cavity modes, which would permit the study of entangle- 
ment between different spin blocks [T?], (iii) controlled 
preparation of robust (insensitive to noise/environment), 
highly entangled states by evolution from an initial 
product state |17l I18j . (iv) measurement of more gen- 
eral atomic spin correlations and their evolution with 
time, which can also provide signatures of criticality in 
QPT's [55 , (v) extending our system to accommodate 
more complex spin models, e.g., by adding additional 
lasers to the setup explained in Sec. |IT] to realize the 
so-called "two-field model" [TS], and (vi) imposing some 
spatial variation on the cavity mode to provide, for ex- 
ample, short ranged interactions, which could be uniform 
or quasirandom. 



FIG. 18: Rescaled concurrence CR(t) for A'^ — 100, with A 
1, r„ = 0.01, and Ft, = 0.2. 



VII. CONCLUSIONS 

We have proposed in this paper a feasible cavity QED 
setup, consisting of a collective atomic pseudospin and 
two quantized cavity modes, which realizes a dissipative 
version of the LMG model in which the interacting spin 
system displays both first- and second-order nonequilib- 
rium quantum phase transitions. The lossy cavity's out- 
put light fields can be utilized to monitor the system as 
the model parameters are varied; specifically, we showed 
that the transmission spectra vary dramatically in the 
vicinity of the transition, with features that are charac- 
teristic of the criticality. A further important result is the 
steady-state entanglement criticality at the QPT and the 
possibility of directly observing this via homodyne detec- 
tion of the cavity output fields. In particular, the entan- 
glement can be quantified rather directly in terms of mea- 
surable atomic quadrature variances. We also observed 
an important sensitivity of the entanglement measure to 
the quadrature phase angle in the critical regimes, which 
we were able to interpret by employing an atomic phase 
space distribution. Finally, we have considered how en- 
tanglement evolves in this system, observing not only 
the criticality at the QPT at long times (corresponding 
to the steady state), but also a rich transient behavior at 
shorter times. 

For future studies, it is clear that the system we 
have proposed offers a variety of opportunities, such as 
(i) investigating phase transitions in response to vari- 
ation of the strength of dissipation (i.e., F;,), (ii) ex- 
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APPENDIX A: COEFFICIENTS OF THE 

ATOM-CAVITY HAMILTONIAN IN THE 

LINEARIZED REGIME 



In Sec. |V A4| we gave the general form of the linearized 
Hamiltonian of the joint atom-cavity system, Eq. (37 1. 



The coefficients of this Hamiltonian in terms of the sys- 
tem parameters, h,Xa,Xb,Tt) and the angles 0,(l> from 
SecEXlare 



Sr. = 



A = 



2h cos 9 + 2 sin 9 [2AXss cos (f> 

-Tb{Yss cos - Xss sin <j))] , 

Xa 



[(1 



+ ( 1 — COS 9) (sin (j) + i cos ( 



Bi 

B2 



Afc 

2 
Afc 
2 



[(1 — cos 0) (sin 0-1- zcos(/))^] 
(1+ cos 61). 



(Al) 

)^]> 
(A2) 

(A3) 
(A4) 



Note that for A < Ac one has = and = 0, giving 
the simplified expressions 5c — 2h, A = Aa, i?i = 0, 
and B2 = Xb- Similar to Sec. |V A2| we can also derive 
simplified expressions in the limit A 3> Ac, i.e., for A — > 
00, one has Sc = 4Aa, A = 0, Bi = —Xb/2, and B2 — Xb/2 
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